function [D, DPol, MargDamagesHour]=calc_damages5(ResProd, Load, RenProd, PPRegion, PPRegion2, PPInt, PPCap, PPInf, PPRandStacked, ...
 ProdMedian,DamagesCoef1,DamagesCoef2, DamagesCoefPoll1,DamagesCoefPoll2) % damages\\

global nhours nnerc npp  PV npol line_losses_loc demand


global PPIntCons PPRegionCons PPRegion2Cons PPInfCons PPRandStackedCons TransCap trans_con_loc

%PPRegion=.001*ones(size(PPRegion));
%PPRegion2=0.*PPRegion2;
%DamagesCoef1=1*ones(size(DamagesCoef1));
%DamagesCoef2=DamagesCoef1;
%PPCap=100000000000+PPCap;

DamagesHighExtra=DamagesCoef2-DamagesCoef1;

ELoad=Load-RenProd-ResProd;


if line_losses_loc==1
    LL_amt=calc_line_losses(demand, ResProd);
    ELoad=demand-RenProd-ResProd+LL_amt;
end


ELoad(ELoad<0)=0;

ELoad2=ELoad.*ELoad;

PPInf2=PPInf.*PPInf;

PPCapSpread=PPCap*ones(1,nhours);

%sum(sum(ResProd))

%sum(sum(ELoad))

%% need to loop over regions to deal with inflection points


NCopies=size(PPRandStacked,3);

DPolStacked=zeros(npol,NCopies);

MargDamagesHourStacked=zeros(nhours, nnerc, NCopies);
          
DStacked=zeros(NCopies,1);

for i=1:NCopies

PPProd=PPInt*ones(1,nhours); % each power plants production by hour over the course of one year.

PPRand=PPRandStacked(:,:,i);

PPProd=PPProd +PPRand;

PPProdMarg=zeros(npp, nhours, nnerc);


if trans_con_loc==1
    PPProdCons=PPIntCons*ones(1,nhours); % each power plants production by hour over the course of one year.
    
    PPRandCons=PPRandStackedCons(:,:,i);
    
    PPProdCons=PPProdCons +PPRandCons;
    
   
    PPInf2Cons=PPInfCons.*PPInfCons;

    [~, MyRegion]=max(PPRegionCons,[], 2);

    MyRegionMat=sparse(MyRegion,1:size(MyRegion,1),1)';
    
    MyRegionMat=full(MyRegionMat);

    MyTransCap=MyRegionMat*TransCap; %npp*1

    MyELoad=MyRegionMat*ELoad'; %npp * nhours

    Exceeds=(MyELoad>(MyTransCap*ones(1,nhours)));

    PPProd(Exceeds)=PPProdCons(Exceeds);


end



for region=1:nnerc
    PPCont=PPRegion(:, region)*ELoad(:, region)' + PPRegion2(:, region)*ELoad2(:, region)'; %npp * nours
    
    PPProdMargReg=PPRegion(:, region)+ 2.*ELoad(:, region)'.* PPRegion2(:, region);
    
   % PPInfCont=(PPRegion(:, region).*PPInf(:, region))*ones(1,nhours) + (PPRegion2(:, region).*PPInf2(:, region))*ones(1,nhours); % prdocution at inflection point. Will
    PPInfCont=(PPRegion(:, region).*PPInf(:, region)) + (PPRegion2(:, region).*PPInf2(:, region)); % prdocution at inflection point. Will
    PPInfCont=PPInfCont*ones(1,nhours);
    %enforce that this is for the correct region
    
   % Select=((ones(npp,1)*ELoad(:, region)')>(PPInf(:, region)*ones(1,nhours))) .* ((PPRegion2(:,region)<0)*ones(1,nhours));
   % Select=Select+((ones(npp,1)*ELoad(:, region)')<(PPInf(:, region)*ones(1,nhours))) .* ((PPRegion2(:,region)>0)*ones(1,nhours));
    Select=((ones(npp,1)*ELoad(:, region)')>(PPInf(:, region))) .* ((PPRegion2(:,region)<0));
    Select=Select+((ones(npp,1)*ELoad(:, region)')<(PPInf(:, region))) .* ((PPRegion2(:,region)>0));   
    
 %   Select=0*Select;
    
    PPProdMargReg(Select==1)=0;
    

    
    PPCont(Select==1)=PPInfCont(Select==1);
    


    if trans_con_loc==1

        PPContCons=PPRegionCons(:, region)*ELoad(:, region)' + PPRegion2Cons(:, region)*ELoad2(:, region)'; %npp * nours
        
        PPProdMargRegCons=PPRegionCons(:, region)+ 2.*ELoad(:, region)'.* PPRegion2Cons(:, region);
        
       % PPInfCont=(PPRegion(:, region).*PPInf(:, region))*ones(1,nhours) + (PPRegion2(:, region).*PPInf2(:, region))*ones(1,nhours); % prdocution at inflection point. Will
        PPInfContCons=(PPRegionCons(:, region).*PPInfCons(:, region)) + (PPRegion2Cons(:, region).*PPInf2Cons(:, region)); % prdocution at inflection point. Will
        PPInfContCons=PPInfContCons*ones(1,nhours);
        %enforce that this is for the correct region
        
       % Select=((ones(npp,1)*ELoad(:, region)')>(PPInf(:, region)*ones(1,nhours))) .* ((PPRegion2(:,region)<0)*ones(1,nhours));
       % Select=Select+((ones(npp,1)*ELoad(:, region)')<(PPInf(:, region)*ones(1,nhours))) .* ((PPRegion2(:,region)>0)*ones(1,nhours));
        SelectCons=((ones(npp,1)*ELoad(:, region)')>(PPInfCons(:, region))) .* ((PPRegion2Cons(:,region)<0));
        SelectCons=SelectCons+((ones(npp,1)*ELoad(:, region)')<(PPInfCons(:, region))) .* ((PPRegion2Cons(:,region)>0));   
        
     %   Select=0*Select;
        
        PPProdMargRegCons(SelectCons==1)=0;
                
        PPContCons(SelectCons==1)=PPInfContCons(SelectCons==1);
        
        PPProdCons=PPProdCons+PPContCons;


        PPProdMargReg(Exceeds)=PPProdMargRegCons(Exceeds);

        PPCont(Exceeds)=PPContCons(Exceeds);


    end

    PPProdMarg(:,:, region)=PPProdMargReg;    

    PPProd=PPProd+PPCont;




end



%dealing with capacity contraints
% have to change marginals first or else these won't turn on
for region=1:nnerc
    PPProdMargReg=PPProdMarg(:,:, region);
    
    PPProdMargReg(PPProd> PPCapSpread)=0;
    
    PPProdMargReg(PPProd<0)=0;
    
    PPProdMarg(:,:, region)=PPProdMargReg;
    
end


PPProd(PPProd<0)=0;

PPProd(PPProd> PPCapSpread)=PPCapSpread(PPProd> PPCapSpread);








%npp*nhours
%DamageMat=PPProd.*(DamagesCoef1 *ones(1,nhours)) + ...
%                  (PPProd-ProdMedian).*(PPProd>ProdMedian).*(DamagesHighExtra *ones(1,nhours)) ;
DamageMat=PPProd.*(DamagesCoef1 ) + ...
                  (PPProd-ProdMedian).*(PPProd>ProdMedian).*(DamagesHighExtra ) ;
    


TotalDamagesYear=sum(DamageMat, 'all');

DStacked(i,1)=TotalDamagesYear*PV;

% power plant by hour by region

          MargPPDamagesHour=zeros(npp, nhours, nnerc);       
for region=1:nnerc              
%    MargPPDamagesHour(:,:, region)=PPProdMarg(:,:, region).*(DamagesCoef1 *ones(1,nhours) + ...
 %                     (PPProd>ProdMedian).*(DamagesHighExtra *ones(1,nhours))) ;
     MargPPDamagesHour(:,:, region)=PPProdMarg(:,:, region).*(DamagesCoef1  + ...
                      (PPProd>ProdMedian).*(DamagesHighExtra)) ;
end

% total damages by region by hour for one year. multiplied by A_t later
% which makes it lifetime. 
MargDamagesHourStacked(:,:,i)=squeeze(sum(MargPPDamagesHour,1));


% Damages by pollutant
for p_i=1:npol
    
    DamagesHighExtraPol=DamagesCoefPoll2(:,p_i)-DamagesCoefPoll1(:,p_i);


    %npp*nhours
%    DamageMatPol=PPProd.*(DamagesCoefPoll1(:,p_i) *ones(1,nhours)) + ...
%                          (PPProd-ProdMedian).*(PPProd>ProdMedian).*(DamagesHighExtraPol *ones(1,nhours)) ;
    DamageMatPol=PPProd.*(DamagesCoefPoll1(:,p_i) ) + ...
                          (PPProd-ProdMedian).*(PPProd>ProdMedian).*(DamagesHighExtraPol ) ;


    TotalDamagesYearPol=sum(DamageMatPol, 'all');

   
    DPolStacked(p_i,i)=TotalDamagesYearPol*PV;
end




end


DPol=mean(DPolStacked,2);
D=mean(DStacked,1);
MargDamagesHour=mean(MargDamagesHourStacked,3);
